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An orbit following code is developed to calculate ion beam trajectories in magnetized plasmas. The 
equation of motion (the Newton's equation) is solved including the Lorentz force term and Coulomb col- 
lisional relaxation term. Furthermore, a new algorithm is introduced by applying perturbation method 
regarding the collision term as a small term. The reduction of computation time is suggested. 



1. Introduction 

An orbit following calculation of charged particles 
is one of the classic problems. It is straightforward, 
but remains to be an important tool for studying par- 
ticle confinement in laboratory plasmas. In high beta 
magnetic confinement devices, Larmor radius of ener- 
getic particles can be comparable to the characteristic 
scale length of the experiment (for example a particles 
in burning plasmas). Guiding center approximation 
fails in the latter cases. pQ [2] 

In this work, the equation of motion is solved 
incorporating Lorentz force term and Coulomb col- 
lisional relaxation term. Since solving the Lorentz 
force term requires much shorter time step compared 
to the guiding center calculation, [3J HI [5J [5] compu- 
tational efficiency is the key. We introduce a new 
algorithm to calculate ion beam trajectories in mag- 
netized plasma by applying perturbation method re- 
garding the Coulomb collisional relaxation term as a 
small perturbation. 

We start our analysis from studying the ion or- 
bital behavior (with and without collision effects) in 
a simple geometry where the magnetic field is axis- 
symmetric. In this paper, we employ a theta pinch 
plasma [71 |H] in a two dimensional system at a plasma 
equilibrium. The orbit following calculation can be 
useful in studying suppression of tilting instabilities [5] 
and rotational instabilities [5] by the ion beams. 

In Sec. m the basic computation model is de- 
scribed. The orbit following calculation is discussed 
in Sec. 121 Sectional presents the perturbation method. 
We summarize this work in Sec. |5l 

2. Equation of motion 

In this section, the equation of motion is de- 
scribed. Ion beam equation in the MKS unit is given 



I— = ffv X B 
dt ^ 



E 



log Aq*' 



dx 
It 



(2) 



Here we recapitulate Ref.^^ as precise as possible 
(including the notations) , for the transparency of the 
work. The first and the second term of Eq.(l) are 
the Lorentz force term (we assume the electric field 
to be zero) and Coulomb collisional relaxation term, 
respectively. The second term reflects the momentum 
change of the test particle per unit time. [TO] Here, 
TO and q are the mass and the charge of the beam 
ions. The magnetic field is given by B while the ion 
beam positions and velocities are given by x and v, 
respectively. The vacuum permittivity is given by Eq, 
and the Coulomb logarithm (see appendix) is given 
by A. All the variables with the superscript * signify 
that of the background plasma species (the ions and 
the electrons). Here, to^ = mm* /{m + to*) is the 
reduced mass. The function $i represents the Gaus- 
sian velocity distribution of the background plasma 
(see appendix), where T* is the background plasma 
temperature and b* = (TO*/2g*r*)^/^. 

Equations (1) and (2) are solved in a Cartesian co- 
ordinate X, y, and z using a fourth order Runge-Kutta- 
Gill method. Equations (1) and (2) holds for ion 
orbital behavior in three dimensional magnetized plas- 
mas in general. In this paper, as an initial application, 
a rigid roter profile of theta pinch plasma [3 [5] is em- 
ployed for the two dimensional magnetic field model. 



Denoting r 

by 



' y2')i/2^ -j-jjg magnetic field is given 
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B = Bo tanh [k {2r'^/rl - l)] z 

where the background density is given by 

n* = nosech^ [n {2r'^/rl - l)] , (4) 

where k is a constant and r^ is the radius at the 
separatrix. [12] Equations (3) and (4) are in plasma 
equilibrium. [7] Correspondingly, the magnetic flux 
-0(r) = B{r')r'dr' is given by 
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log [cosh [K(2rVr^ - 1)]] 
log [cosh (k)]. (5) 



In this paper, the angular momentum is given bv|10) 



Pq = mr^O + q'il!{r), 



(6) 



where 9 is the time derivative of the angular coordi- 
nate 9. The kinetic energy is given by Ek — mv^/2. 

3. Beam ion orbit 

In this section, the ion orbit calculation is pre- 
sented employing Eqs.(l) and (2). We study beam 
ion (energetic particle) behavior whose temperature 
is much larger than that of the background thermal 
plasma. 

Figure 1 shows the particle orbits in a Cartesian 
coordinate in the absence of Coulomb collisions. The 
magnetic configuration reflects that of the FRC injec- 
tion experiment (FIX) parameter in the confinement 
chamber; [13] in Eq.(3), the magnetic field strength Bo 
is 0.05(T) and the separatrix radius is given by 
0.2(m) and thus the magnetic null is at r„ — 0.141(m). 
The wall radius is set at r^, = 0.4(m). In Eq.(3), 
we set K = 0.6. [7] The beam ion species is Hydro- 
gen. Throughout this paper, we assume that the 
neutral beams are ionized at x = — 0.136(to) and 
y = 0.147(m) which is on the separatrix (the initial 
position of the beam ion calculation is given there). 

In Fig. 1(a), the beam ion temperature is given 
by Tb = 50{eV) [followed the ion orbit for 20(^^5)], 
while in Fig. 1(b), the beam ion temperature is given 
by Tb = 4000(ey) [followed the ion orbit for 2(/zs)]. 
Naturally, the Fig. 1(b) case has a larger Lamor ra- 
dius. As one can see, the direction of the Larmor pre- 
cession changes when the trajectory crosses the mag- 
netic null point "r„" . This is referred to as meandering 
motion. ' Si Since the magnitude of the magnetic field 
inside r„ is weaker than the outside, the Larmor ra- 
dius is slightly larger inside the separatrix [see Fig. 1(c) 
where the magnetic field strength and the density pro- 
file are depicted] . As shown above the motion is peri- 
odic which can be understood by the Nocther's theo- 
rem (a canonical variable conjugate to a constant mo- 
mentum undergoes periodic motion). The conserva- 
tion of kinetic energy and the angular momentum is 
verified for the calculation in Fig.l. With a single 
precision, the momentum (energy) conserves at the 
accuracy of 6.3 x 10-^% (5.5 x 10"^%) of the absolute 
value after following the orbit for 10000 steps. Here, 
the time step in the calculation is given by one per- 
cent of 2tt/Qc where Qc is the beam ion's cyclotron 
frequency. 

Figures 2 and 3 show the particle orbits in the 
presence of Coulomb collisions. The collision effect 
is dominated by electrons (see appendix). The back- 
ground electron density and temperature is given by 
no = 5 X 10^^{m'^) and = 20{eV). In Fig.2, the 
beam ion temperature is given by Tb = lOO(ey). In 
Fig.3, Tb = 2000(ey). 
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Fig. 1 Orbital behavior of (a) 50 (eV) (b) 4000 (eV) Hy- 
drogen ions. The collision term is turned off in 
Eq.(l). The Larmor precession changes its direc- 
tion at the magnetic null point (dashed circle) . Unit 
length is normalized by Vs (solid circle), (c) The 
magnetic field strength (green curve) and the den- 
sity profile (red curve) are depicted. The location 
of the separatrix and the magnetic null point are 
suggested. 




2 (a) Orbital behavior of 100 eV beam ion in the pres- 
ence of the coUision term. The background electron 
density and temperature is given by 5 X 10^^ (m-^) 
and 20{eV). (b) The kinetic energy (black curve) 
and the angular momentum (red curve) versus 
time. 



Figures 2(b) and 3(b) show the kinetic energy 
(Ek) and the canonical angular momentum {P0) ver- 
sus time. In Fig. 2 and 3, the kinetic energy and the 
angular momentum relax. The e-fold times estimated 
in Fig. 2(b) and Fig. 3(b) are summarized in Table 1 
(we take the logarithm; r*."" for the kinetic energy 
and T^"^ for the angular momentum). 

We now compare the numerical relaxation time [e- 
fold time estimated from Fig. 2(b) and Fig. 3(b)] with 
a theoretically estimated relaxation time. [TUJ [T3] Fol- 
lowing Ref . [lOj , the energy relaxation time is given by 

Tie log Aq^'q^ \m rriej 



Table 1 : Comparison of relaxation times. 
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(which is independent of beam ion temperature unless 
Tb/m ~ Te/nie) and the perpendicular momentum 
relaxation time is given by 



nel0gAg2gr2$(5^^)' 

respectively. The background plasma is assumed to 

be only electrons. In Eqs.(7) and (8), the charge, the 
mass, and the temperature of electrons are given by qe, 
irie, and Tg, respectively. Here, be = (me/2geTe)^^^. 
The relaxation time employing the parameters used 
in Fig. 2 and Fig. 3 are summarized in Table 1. In Ta- 
bid, the energy relaxation time compares favorably 
with the numerical estimation, while the momentum 
relaxation time differs in particular for the higher en- 
ergy case. 

4. Perturbation method 

In this section, the perturbation method is intro- 
duced. Normalizing Eqs.(l) and (2) by the beam ion 
cyclotron frequency Clc = qbBo/nib, and the separatirx 
radius Tg, we obtain 



— = V X B - eF 

dT 



dK 

dT 



V 



(9) 



(10) 



where the frictional force is regarded as a small term 
employing 



gMogAno $i(&*) 

e = - — 5 ^TTj > <C 1 
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and 



Expanding B = Bq + eBi for the rigid rotor profile, 
we have 

Bo = tanh [k {2R^ - l)] z (13) 

Bi = 4k(Xo • Xi)sec/i? [n {2R^ - l)] z (14) 

N{R) = sech"^ [k {2B? - 1)] (15) 

Here, the capital letters (T, V, i?, and X) represent 
the normalized time, velocity, radius, and position, 
respectively. 



From Eqs.(9) and (10), the lowest order equation 
is given by 



dT 
dXo 

dT 



Vo X Bo 



Vo 



and the first order equation in order e is given by 
dVi 



dT 
rfXi 

It 



= Vi X Bo + Vo X Bi -h F 
= Vi. 



(16) 
(17) 

(18) 

(19) 



The solution then is given by the summation X = 
Xo -I- eXi, V = Vo eVi. The crux in Eqs.(18) and 
(19) arc the changes in particle velocity (Vi x Bo) and 
particle's displacement (Vo x Bi) both induced by the 
small friction force (the F term). 

The perturbation method is useful since we only 
need to change the constant e when the the plasma pa- 
rameters change, e.g. background densities and tem- 
peratures (and do not need to recalculate the whole 
trajectories). Figure 4 and 5 show particle trajecto- 
ries when the perturbation method is employed. Here, 
the green curve solution in Fig. 5 are obtained by re- 
cycling Xo and Xi from Fig. 4, by simply changing 
the parameter e. In both Figs. 4 and 5, the solution 
from the perturbation method [green curves, solved 
Eqs.(16)-(19)] matches with the direct collision calcu- 
lations [red cm-vcs, solved Eqs.(l) and (2)]. 

The lowest order solution is periodic when the 
magnetic field is axis-symmetric. Likewise we expect 
the first order solution to be periodic. If the lattcir is 
the case, there will be another attractive application 
of the perturbation method. By storing the first pe- 
riodic motion of both the lowest and the higher order 
solution, the algorithm can predict periodic motion 
in the later phase and thus can reduce computation 
time. As a demonstration, here we take a simplified 
case where the Lorentz force and the centrifugal force 
are balanced at the initial state; 



mrO^ = qvB 



(20) 



[the trajectory will be a perfect circle in the absence 
of collisions. See Fig. 6(a)]. Figure 6(b) suggests a 
periodic motion of the first order solution from the 
perturbation method (time evolution of the Cartesian 
coordinates x\ and y\ are plotted). 

5. Summary 

An orbit following code is developed to calculate 
ion beam trajectories in magnetized plasmas. The 
equation of motion is solved incorporating the Lorentz 
force term and Coulomb collisional relaxation term. 
Conservation of energy and angular momentum is con- 
firmed in the absence of collisions. With the collisions. 
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(a) The beam ion orbit with a background plasma 
Te = 50eV and = 1.0 x 10^^. The beam ion 
temperature is 300(6^). (b) Expansion of the final 
stage of Fig. 4(a). The green (red) curve are from 
the perturbation method (the direct calculation). 



Fig. 5 (a) The beam ion orbit with a background plasma 
Te = 50e1/ and = 5.0 x 10^^. The beam ion 
temperature is 300(el/). (b) Expansion of the final 
stage of Fig. 5(a). Here, the green curve solution in 
(b) are obtained by recycling Xq and X\ from (a), 
by simply changing the parameter e. The green 
(red) curve are from the perturbation method (the 
direct calculation). 
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Fig. 6 (a) The beam ion orbit (red circle) when the 
Lorentz force and the centrifugal force are balanced 
at the initial state. The initial position is at a; = 
and y = Vs. (b) Time evolution of the perturbed 
quantities xi (solid), j/i (dashed). 



it is shown that the energy relaxation time compares 
favorably with the theoretical prediction. [TU] 

Furthermore, a new algorithm to calculate ion 
beam trajectories is reported. We have applied per- 
turbation method regarding the collisional term small 
compared to the zeroth order Lorentz force term. The 
two numerical solutions from the perturbation method 
and the direct collisional calculation matched. The 
perturbation method is useful since we only need to 
change the perturbation parameter to recalculate the 
trajectories, when the background parameters change. 
In general, the algorithm can be applied to periodic 
motion under perturbative frictional forces, such as 
guiding center trajectories |31 HI [5J [B] or satellite mo- 
tion. We have also suggested a reduction in compu- 
tation time by capturing the periodic motion. More 
detailed analysis will be our future work. 

The work is initiated as a diploma thesis at 
Osaka University during the years 1990-1991. [TS| The 
author would like to thank Dr. T. Ishimura and Dr. 
S. Okada for useful discussions. A part of this work 
is supported by National Cheng Kung University Top 
University Project. The author would like to thank 
Dr. C. Z. Cheng and Dr.K. C. Shaing. 

Appendix 

We recapitulate Coulomb collision processes pre- 
sented in Ref.[Tn]. Assuming the background plasma 
is Maxwellian, change in the momentum for the beam 
ion within time dt is given by 



dp 

It 
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logA(g*) 



■n*$i(6*w)(21) 



here, q (q*) and m (m*) are the charge and the 
mass of the beam ions (background plasma species), 
rrir = mm* /{m + m*) is the reduced mass, n* is the 
background plasma density. Here, ^* signifies sum- 
mation over species. The Coulomb logarithm is given 
by 



log A ~ 7 -I- log 
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Tie \ 
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1/2 



(22) 



where the electron temperature is in the unit of eV , 
and the electron density is in the unit of to^"^. 

Letting X ^ h*v ^ {m* /2eT*f''^ {2eT /mf'^ = 

1/2 

{m*T/mT*) ' , the function ^{x) is given by 
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$i(a;) = $(x) j= exp i 



(23) 
(24) 



In Eq.(21), the contribution from electrons is much 
larger than the ions because of the 1/mr factor. The 
relaxation of high energy ions is namely due to the 



collision with the electrons. One can then employ the 
form 

= ^. (25) 

at the X ^ 1 limit [the integration form of Eqs.(23) 
and (24) are employed in the computation]. 



References 

[1] E. R. Solano, Phys. Plasmas 3, 1187 (1996). 

[2] D. R. Mikkelsen, Phys. Plasmas 4, 3667 (1997). 

[3] A. I. Morozov and L. S. Solevev, Reviews of Plasma 

Physics Vol.2, (Consultants Bureau, New York, 

1966), P.201. 

[4] R. G. Littlcjohn, J. Plasma Phys. 29, 111 (1983). 
[5] Y. Nishimura and M.Azumi, Phys. Plasma 4, 2365 
(1997). 

[6] Y. Nishimura, Contributions to Plasma Physics 48, 
224 (2008). 

[7] R. L. Morse and J. P. Preidberg, Phys. Fluids 13, 531 
(1970). 

[8] T. Ishimura, Phys. Fluids 27, 2139 (1984). 

[9] R. Horiuchi and T. Sato, Phys. Fluids B 2, 2652 

(1990) . 

[10] K. Miyamoto, Plasma Physics for Nuclear Fusion 2nd 

ed., (Iwanami, Tokyo, 1986), p. 76. 
[11] I. Kawakami, Suuchi-keisan, (Iwanami, Tokyo, 1989), 

p. 159 (in Japanese). 
[12] M. Tuszewski, Nucl. Fusion 28, 2033 (1988). 
[13] A. Shiokawa and S. Goto, Phys. Fluids B 5, 534 

(1993). 

[14] L. Spitzer, Jr., Physics of Fully Ionized Gases, (In- 

terscience. New York, 1962). 
[15] Y. Nishimura, Diploma thesis, Osaka University 

(1991) . 



